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Considerations  for  Stationary 
Ice  Covered  Flows 
in  Adaptive  Hydraulics  (ADH) 


by  Gary  L.  Brown,  Gaurav  Savant,  Charlie  Berger,  and  David  S.  Smith 


PURPOSE:  This  System-Wide  Water  Resources  Program  (SWWRP)  technical  note  presents 
the  theoretical  development  and  implementation  of  the  equations  and  techniques  necessary  to 
simulate  stationary  ice-covered  flow  in  ADaptive  Hydraulics  (ADH).  This  implementation 
includes  the  application  of  a  surface  pressure  field  to  simulate  the  weight  of  floating  ice  cover  on 
the  flow,  as  well  as  a  method  developed  at  the  U.S.  Army  Engineer  Research  and  Development 
Center  (ERDC)  for  the  calculation  of  drag  resulting  from  combined  bed  roughness  and  ice  cover. 

BACKGROUND:  The  presence  of  ice  cover  over  flowing  water  complicates  the  hydraulic 
properties  of  the  flow  by  modifying  the  available  flow  cross-sectional  area  (via  the  floating  ice 
cover),  and  by  modifying  the  resistance  to  flow  (by  increasing  the  surface  area  available  to  shear 
forces  and  contributing  unique  roughness  characteristics  associated  with  the  ice  cover  to  the 
development  of  the  shear  profile)  (Ashton  1986).  The  two-dimensional  shallow-water  (SW2) 
module  of  ADH  has  been  extended  to  include  changes  to  the  flow  due  to  the  pressure  imposed 
by  an  ice  field  as  well  as  modifications  to  the  shear  forces.  The  techniques  employed  to  account 
for  these  changes  are  described  in  detail  in  this  technical  note. 

SHEAR  STRESS  INDUCED  BY  BOTH  BED  AND  ICE  COVER:  In  general,  the  presence 
of  ice  cover  significantly  affects  the  hydraulic  roughness  associated  with  a  flow.  If  one  applies 
classic  hydraulic  velocity  profile  theory  to  the  ice  cover  problem,  the  complete  velocity  profile 
can  be  represented  by  two  hydraulically  independent  profiles,  which  share  a  single,  maximum 
velocity  (for  a  review  of  classical  hydraulic  theory,  see  Schlichting  1979).  These  profiles  divide 
the  total  flow  into  two  distinct  regions,  the  ice  region  and  the  bed  region.  The  ice  region  and  the 
bed  region  are  the  flow  regions  above  and  below  the  maximum  velocity  line,  respectively 
(Figure  1). 

Several  researchers  have  questioned  the  validity  of  this  two-layer  theory  (e.g.,  Ettema  2002).  The 
classical  theory  dictates  that  the  eddy  viscosity  coinciding  with  the  maximum  velocity  should  be 
exactly  zero,  but  observations  indicate  this  eddy  viscosity  is  generally  non-zero  (Krishnappan 
1983).  Also,  observations  of  the  shear  stress  profile  indicate  that  the  location  of  zero  shear  stress 
does  not  necessarily  coincide  with  the  maximum  velocity  (Parthasarathy  and  Muste  1994).  These 
observations  imply  that,  contrary  to  the  claims  of  classical  theory,  some  energy  is  exchanged 
between  profiles,  and  hence,  the  profiles  are  not  hydraulically  independent. 
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Figure  1 .  Schematic  Representation  of  the  Velocity  and  Shear  Profiles  under  Ice  Cover 


In  order  to  develop  a  new  expression  for  shear  stress  under  ice  cover,  the  classical  hydraulic 
theory  is  assumed  valid  and  the  classical  theory  for  turbulent  rough  flow  between  two  parallel 
surfaces  with  different  hydraulic  roughnesses  is  developed.  An  approximation  term  is  added  for 
the  cross-profile  exchange  of  turbulent  shear  by  assuming  that  this  quantity  is  small  relative  to 
the  total  drag  and  can  be  approximated  with  a  linear  superposition  of  a  correction  tenn  derived 
from  the  classical  theory. 

Approach.  First,  the  classical  theory  for  turbulent  rough  flow  between  two  parallel  surfaces 
with  different  hydraulic  roughnesses  is  developed.  The  following  three  assumptions  arise  from 
the  classical  theory: 

1.  The  ice  region  and  bed  region  have  the  same  maximum  velocity,  which  is  located  at  the 
junction  of  the  profiles. 

2.  Since  the  vertical  velocity  gradient  at  the  maximum  velocity  height  is  zero,  the  shear  stress  at 
the  maximum  velocity  height  must  also  be  zero. 

3.  Since  no  energy  can  be  exchanged  across  a  horizontal  plane  of  zero  shear  stress,  the  two 
velocity  profiles  are  hydraulically  independent.  Therefore,  the  total  shear  is  equal  to  the  sum 
of  ice  and  bed  shear  stresses,  which  can  be  evaluated  independently. 

Tx  ~  TBED.x  tice.x  (1) 


Ty  ^BED.y  ^ICE.y 


The  shear  stresses  in  x-  and  y-directions  are  given  as  follows: 


TBED.x  —  2  PRO  BED (aBEDVx  )'\/(aBEDVx  )  (aBEDVy) 


ttBEDVy  /  "I"  TCPE.x 


TBED.y  —  2  P^D  BED  (ttBED  Vy  )\/(aBEDVx)  "*"(aBEDVy)  T 


CPE.y 


(2) 

(3) 

(4) 
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TICE.x  2  P^D.ICE(aiCEVx)'\/(aiCEVx)”  "*~(aiCEVy)  T 
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TICE.y  —  2  P^D.ICE  (aiCEVy  \j (ttICEVx  )  (aiCEVy  )  1 

~~  ^(PbedVi 
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:-ln(piCE)uf 
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MCE 


f.ICE 


(10) 


PBED  =29.7— ^  +  1 


‘'BED 


P,ce=29.7^=J  +  1 

k, 


"ICE 


k  =  0.4 


(ID 

(12) 

(13) 


Variables  are  defined  at  the  end  of  this  technical  note. 

Note  that  Equations  11  and  12  are  expressions  of  the  form  of  the  velocity  profile  given  by 
Christensen  (1972).  This  velocity  profile  is  identical  to  the  traditional  velocity  profile,  except  for 
the  addition  of  the  “+1”  term.  This  additional  term  only  has  a  significant  effect  on  the  profile  for 
small  values  of  the  roughness  ratio  (z/k).  It  ensures  that  the  velocity  magnitude  v  is  greater  than 
zero  for  all  possible  values  of  z/k  (the  traditional  profile  approaches  -oo  for  z/k  — >  0). 

According  to  the  classical  theory,  the  shear  profile  is  linear  for  both  the  bed  and  ice  profiles,  and 
the  shear  at  zmv  must  equal  zero  (since  the  velocity  at  zmv  is  the  maximum  velocity  and  hence,  the 
inflection  point  of  the  total  velocity  profile).  This  assumption  can  be  used  to  express  the  ratio  of 
the  bed  and  ice  shear  stresses  as  a  linear  function  of  zmv  and  d. 
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Consider  the  force  balance  for  the  total  depth: 

*  =  PgdSEGL 

Now,  consider  the  force  balance  for  just  the  bed  friction: 


TBED  P§Zmv^EGL  (  1 5) 

Equation  15  must  be  the  equation  for  the  bed  force,  because  no  energy  transfer  can  occur  across 
the  zero  shear  stress  line,  which  occurs  at  zmv.  Hence,  the  bed  shear  can  only  resist  the  weight  of 
the  water  from  zero  to  zmv.  The  remainder  of  the  energy  (from  zmv  to  d)  is  absorbed  by  the  ice 
shear.  Next,  note  that: 


T  —  T  bed  T  ICE 

Substituting  Equation  16  into  Equation  14  and  divide  Equation  14  by  Equation  15  yields: 


(16) 


Zmv  _  TBED 
d  —  Zmv  ^  ICF 

Combining  Equations  7,  8,  and  17  can  obtain  an  equation  for  zmv. 


(17) 


Zmv  {^7T)d 

where: 

_  ^(Pbed) 

ln(pICE) 


(18) 


(19) 


Equations  18  and  19  must  be  solved  iteratively,  but  the  solution  converges  rapidly.  A  good  initial 
estimate  of  zmv  to  begin  the  iterations  is  zmv  =  d/2. 

The  location  of  zmv  is  now  known  and  is  used  to  determine  the  ice  and  bed  drag  coefficients. 
They  are  found  by  solving  Equations  3-6  for  the  drag  coefficients,  invoking  Equations  9  and  10 
to  express  the  shear  stresses  in  terms  of  friction  velocities,  and  then  integrating  the  velocity 
profile  to  derive  an  expression  for  av/uf  in  tenns  of  p. 


C  =2 

^D.BED 


K 


(Pbed-1) 


vtPBED{ln(pBED)-1}+1] 


(20) 
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C 


D.ICE 


K(Picn  —  0 

(j3icE{ln(PicE)-li+l]J 


(21) 


The  mean  velocity  correction  factors  adjust  for  the  ratio  between  the  mean  velocities  of  each 
profile  individually  and  the  mean  velocity  over  the  entire  water  depth.  These  relationships  are 
derived  from  Equations  3-6,  Equations  9  and  10,  as  well  as  the  ratio  given  in  Equation  17. 
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r 
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Now  that  the  classical  theory  has  been  developed,  an  approximate  method  can  be  derived  from 
the  classical  theory  to  account  for  the  energy  transfer  between  profiles.  The  theoretical 
development  given  here  is  derived  with  the  implicit  assumption  that  no  energy  can  be  exchanged 
between  velocity  profiles.  However,  several  researchers  have  noted  that  the  eddy  viscosity  is 
generally  nonzero  at  the  profile  interface  (Krishnappan  1983).  Figure  2  is  a  schematic  of  the 
difference  between  the  theoretical  and  observed  eddy  viscosity  profiles. 


Figure  2.  Schematic  Representation  of  the  Eddy  Viscosity  Profile  under  Ice  Cover. 


The  nonzero  eddy  viscosity  at  the  profile  interface  can  cause  some  energy  exchange  across  the 
interface.  Further,  the  differences  in  the  bulk  momentum  of  each  profile  will  cause  some 
momentum  transfer  via  the  finite  mixing  length  associated  with  the  nonzero  eddy  viscosity  at  the 


5 


ERDC  TN-SWWRP-09-4 
May  2009 


profile  interface.  This  contribution  is  exactly  zero  when  both  profiles  are  identical,  and  increases 
as  the  difference  in  the  profiles  increases.  If  one  makes  the  simplifying  assumption  that  this  bulk 
momentum  transfer  across  profiles  can  be  superimposed  onto  the  existing  theory  (which  was 
developed  assuming  no  energy  transference),  the  cross-profile  energy  exchange  can  be  accounted 
for  in  an  approximate  sense  with  the  following  relations. 


CPE.x 


=  P 


SMAX.BED  eMAX.ICE 

(aBED  ^ice  ) 

l  2  J 

l  d(l-5j  J 

1  CPE.y 
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1 
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^  MAX. ICE 


1 

=  —  KU 
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f.ICE 


(d-zmv) 


(28) 


5mv  =  0.368 


(29) 


Note  that  this  cross-profile  exchange  of  momentum  has  no  effect  on  the  total  combined  bed  and 
ice  shear  stress,  since  it  only  passes  momentum  from  one  profile  to  the  other.  Therefore,  the  total 
shear  is  fully-described  by  the  classical  development,  and  only  the  partitioning  of  the  shear 
between  the  bed  and  the  ice  cover  is  affected. 

Comparison  to  Observed  Data.  Flume  experiments  for  simulated  ice  cover  were  conducted 
by  Parthasarathy  and  Muste  (1994).  Calculations  using  the  proposed  ERDC-Coastal  and 
Hydraulics  Laboratory  (CHL)  method  were  compared  to  these  experimental  results  for  several 
values  of  the  ice-to-bed  roughness  ratio.  The  results  are  given  in  Figure  3. 

The  comparisons  show  that  the  proposed  method  is  in  good  agreement  with  the  experimental 
values  for  the  total  shear.  The  comparison  for  the  partition  of  the  shear  stress  between  the  bed 
and  ice  shear  is  not  as  good,  which  may  indicate  that  the  approximation  of  the  cross-profile 
exchange  could  be  improved.  However,  the  experimenters  noted  that,  for  the  case  with  the 
highest  ice-to-bed-roughness  ratio,  they  observed  significant  secondary  currents.  They  postulated 
that  these  currents  may  have  arisen  from  sidewall  effects  in  the  flume  experiments.  Such  currents 
would  increase  the  cross-profile  exchange  of  momentum.  Hence,  the  difference  between  the 
predicted  and  observed  values  may  be  partially  attributable  to  the  additional  mixing  introduced 
by  the  secondary  currents.  Further  experimentation  and  analyses  are  needed  to  quantify  the 
accuracy  of  the  cross-profile  exchange  approximation. 
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Comparsion  Between  Proposed  ERDC-CHL  Method  and  Experimental  Results 
(Parthasarathy  and  Muste,  1994) 


X  Cd  bed  ERDC-CHL 
X  Cd  ice  ERDC-CHL 
+  Cd  total  ERDC-CHL 
□  Cd  bed  Experiment 
o  Cd  ice  Experiment 
A  Cd  total  Experiment 


k  ice  /  k  bed 


Figure  3.  Comparison  of  the  Proposed  ERDC-CHL  Method  to  the  Experimental  Data  of  Parthasarathy 
and  Muste  (1994) 


Expression  for  Manning’s  n.  Traditionally,  the  shear  stress  has  often  been  expressed  in 
terms  of  Manning’s  n.  For  example,  the  Sabaneev  equation  (Nezhikovskiy  1964)  gives  an 
expression  for  Manning’s  n  for  the  composite  ice-bed  roughness  as  a  function  of  Manning’s  n  for 
the  bed  roughness  alone.  The  equation  is  derived  by  invoking  both  velocity  profile  relationships 
and  Manning’s  equation.  It  is  effectively  a  simplified  form  of  the  more  robust  development  given 
by  Larsen  (1969).  The  Sabaneev  equation  is  given  as  follows: 


n  COMPOSITE 
nBED 


V/2 

11  BED 

2n 


+  n 

3/2 


3/2  A 
ICE 


(30) 


An  analogous  relationship  for  the  method  developed  in  this  technical  note  is  given  as  follows: 


n  COMPOSITE 
nBED 


a2  +  C  a2 

V-'D.BEDUBED  ^  V-'  D.ICEUICE 


2C 


D.BED.ONLY 


(31) 


where: 

Cd.bed.only  is  Cd.bed  calculated  for  zmv  =  d/2. 
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IMPOSITION  OF  PRESSURE  FIELD  TO  SIMULATE  FLOATING  ICE  COVER:  The 

presence  of  ice  cover  imposes  pressure  in  addition  to  the  fluid  hydrostatic  pressure  on  the  fluid 
underneath.  This  pressure  is  a  function  of  the  ice  thickness,  as  well  as  the  ice  density 
(Equation  32). 

Pice  =  Pice  Shoe  (32) 

IMPLEMENTATION  INTO  ADH:  Implementing  stationary  ice  effects  on  fluid  hydraulics  in 
ADH  was  done  by  specifying  a  card  in  the  ADH  boundary  conditions  file  (*.bc)  using  one  of  two 
approaches  -  by  material  type  or  by  a  radius. 


To  specify  the  presence  of  ice  by  material  type,  one  must  use  three  cards,  FR  ICE,  FR  IRH,  and 
FR  BRH  (Tables  1,  2,  and  3,  respectively). 


Table  1 

FR  ICE  Card  Description 

Field 

Type 

Value 

Description 

1 

Character 

FR 

Card  type 

2 

Character 

Ice 

Card  type 

3 

Integer 

>  0 

Ice  material  string 

4 

Real 

# 

Ice  thickness 

5 

Real 

# 

Ice  density 

6 

Integer 

0 

Stationary  ice 

Table  2 

FR  IRH  Card  Description 

Field 

Type 

Value 

Description 

1 

Character 

FR 

Card  type 

2 

Character 

IRH 

Card  type 

3 

Real 

# 

Ice  roughness  height 

Table  3 

FR  BRH  Card  Description 

Field 

Type 

Value 

Description 

1 

Character 

FR 

Card  type 

2 

Character 

BRH 

Card  type 

3 

Real 

# 

Bed  roughness  height 
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To  specify  the  presence  of  ice  by  overlapping  circular  regions,  the  INS  (ice  node  string)  card 
must  be  included  (Table  4),  in  addition  to  the  three  cards  used  when  specifying  ice  by  material 
type. 


Table  4 

INS  Card  Description 

Field 

Type 

Value 

Description 

1 

Character 

INS 

Card  type 

2 

Real 

# 

X  coordinate  of  circle  center 

3 

Real 

# 

Y  coordinate  of  circle  center 

4 

Real 

# 

Radius  of  ice  circle 

5 

Integer 

0 

Stationary  ice 

ADDITIONAL  INFORMATION:  For  additional  information,  contact  Gary  L.  Brown,  Coastal 
and  Hydraulics  Laboratory,  U.S.  Army  Engineer  Research  and  Development  Center,  3909  Halls 
Ferry  Road,  Vicksburg,  MS  39180  at  601-634-3628,  e-mail:  Gary.L.Brown (Susace.army.mil,  or 
Dr.  Gaurav  Savant,  P.E.,  Coastal  and  Hydraulics  Laboratory,  U.S.  Army  Engineer  Research  and 
Development  Center,  3909  Halls  Ferry  Road,  Vicksburg,  MS  39180  at  601-634-3628,  e-mail: 
Gaurav.Savant@usace.army.mil.  This  effort  was  funded  through  the  System-Wide  Water 
Resources  Program  (SWWRP).  For  information  on  SWWRP,  please  consult 
https://swwrp.usace.army.mil/  or  contact  the  Program  Manager,  Dr.  Steven  L.  Ashby: 
Steven.L.Ashbv@usace.army.mil.  This  SWWRP  technical  note  should  be  cited  as  follows: 

Brown,  G.  L.,  G.  Savant,  C.;  Berger,  and  D.  S.  Smith.  2009.  Considerations  for 
stationary  ice  covered  flows  in  ADaptive  Hydraulics  (ADH)  ERDC  TN-SWWRP- 
09-4.  Vicksburg,  MS:  U.S.  Army  Engineer  Research  and  Development  Center. 

https  ://swwrp.  usace.  army.  mil/. 
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VARIABLES 


Cd.bed  = 
Cd. bed.  only  = 
Cd.ice 
d  = 

g  = 
kBED  = 
klCE  - 

tlBED  = 
n  COMPOSITE  = 
nicE 

p  = 

1  ice 

lice 

V  = 
V.r  = 
Vv  = 
VMAX.BED  = 

Vmax.ice  = 

Ilf  BED  = 
UfICE  = 
Zmv 


O.BED  ~ 

mice  = 

O-IBR  = 

SMAX.BED  = 
SMAX.ICE  = 

3mv 

K  = 

P  = 

Pice  — 

Segl  = 

TBED  = 
Tice  = 
tcpe 

;  = 


the  bed  shear  stress  drag  coefficient 

the  bed  shear  stress  drag  coefficient,  assuming  zmv  =  d/2 

the  ice  shear  stress  drag  coefficient 

the  water  depth 

the  gravitational  acceleration 

the  equivalent  bed  roughness  height 

the  equivalent  ice  roughness  height 

the  Manning’s  n  for  the  bed 

the  Manning’s  n  for  the  combined  ice-bed  roughness 
the  Manning’s  n  for  the  ice 
pressure  induced  by  the  ice  cover 
ice  thickness 

the  depth-averaged  velocity  magnitude 
the  depth-averaged  velocity  in  the  x-direction 
the  depth-averaged  velocity  in  the  y-direction 
the  maximum  velocity  for  the  bed  profile 
the  maximum  velocity  for  the  ice  profile 
the  friction  velocity  for  the  bed  profile 
the  friction  velocity  for  the  ice  profile 

the  depth  at  which  the  maximum  velocity  is  located  (i.e.  the  location  of  the 
transition  from  the  bed  induced  velocity  profile  to  the  ice  induced  velocity 
profile) 

the  mean  velocity  correction  factor  for  the  bed  shear  stress 

the  mean  velocity  correction  factor  for  the  ice  shear  stress 

the  ratio  of  the  mean  velocity  for  the  ice-induced  velocity  profile  to  the  mean 

velocity  for  the  bed-induced  velocity  profile 

the  maximum  eddy  viscosity  for  the  bed  profile 

the  maximum  eddy  viscosity  for  the  ice  profile 

the  normalized  fraction  of  the  distance  to  the  centroid  of  the  velocity  profile 

the  Von  Kannan  constant 

the  density  of  water 

the  density  of  ice 

the  slope  of  the  energy  grade  line 

the  boundary  shear  at  the  flow-bed  interface 

the  boundary  shear  at  the  flow-ice  interface 

the  approximate  cross-profile  exchange  of  shear  stress 

the  total  shear  stress 


NOTE:  The  contents  of  this  technical  note  are  not  to  be  used  for  advertising,  publication,  or 
promotional  purposes.  Citation  of  trade  names  does  not  constitute  an  official  endorsement  or 

approval  of  the  use  of  such  products. 
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